/****************************************************************\
*                                                                *
*  Codon Substutition Matrix Object                              *
*                                                                *
*  Guy St.C. Slater..   mailto:guy@ebi.ac.uk                     *
*  Copyright (C) 2000-2009.  All Rights Reserved.                *
*                                                                *
*  This source code is distributed under the terms of the        *
*  GNU General Public License, version 3. See the file COPYING   *
*  or http://www.gnu.org/licenses/gpl.txt for details            *
*                                                                *
*  If you use this code, please keep this notice intact.         *
*                                                                *
\****************************************************************/

/* Codon Substitution Data taken from:
 * "Empirical codon substitution matrix"
 * Adrian Schneider, Gina M Cannarozzi and Gaston H Gonnet
 * BMC Bioinformatics. 2005, 6:134
 */

/* In alphabetical order of codon */

#include <stdio.h>
#include <string.h>
#include <ctype.h> /* For tolower() */
#include <math.h>  /* For log()     */

#include "codonsubmat.h"

static gint local_codon_substitution_count
                  [CODON_ALPHABETSIZE_ACGT]
                  [CODON_ALPHABETSIZE_ACGT] = {
/* AAA */
{192416,   5180, 139387,   5795,   4300,   1732,   1129,   1860,
  20646,   2819,  12905,   2447,   1031,    476,   1993,    580,
   7019,   1702,  10180,   1597,   1178,    530,    358,    689,
   4714,   4108,   5511,   2693,    363,    366,    884,    398,
   8048,   1636,   5765,   1836,   1979,   1019,    494,   1188,
   1686,    799,    938,    654,    591,    364,    976,    468,
      0,    247,      0,    261,   1159,    733,    292,    981,
      0,    280,    264,    298,    341,    143,    442,    172},
/* AAC */
{  5180, 142704,   5105,  76958,   2579,   4969,    943,   3302,
   1338,  19096,   1159,   9374,    309,   1115,    772,    663,
   1484,   5608,   3282,   3172,    567,    801,    232,    663,
    333,    942,    541,    433,    126,    446,    471,    263,
   2525,  10085,   3034,   6768,   1131,   2290,    489,   1536,
   1728,   4186,   1544,   1992,    242,    636,    590,    466,
      0,   1308,      0,    803,   1077,   2021,    360,   1682,
      0,    817,    100,    575,    139,    395,    238,    305},
/* AAG */
{139387,   5105, 246420,   5132,   3195,   1797,   1538,   1577,
  15336,   2826,  19246,   2259,    649,    503,   3109,    515,
   5156,   1891,  15687,   1526,    945,    524,    430,    604,
   4333,   5231,   7539,   2968,    281,    370,   1292,    364,
   5401,   1605,   9352,   1732,   1492,   1138,    599,   1064,
   1215,    885,   1386,    628,    432,    370,   1230,    389,
      0,    273,      0,    257,    919,    702,    298,    830,
      0,    311,    441,    290,    263,    120,    530,    170},
/* AAT */
{  5795,  76958,   5132, 107878,   2735,   2939,    939,   4546,
   1376,  11367,   1093,  13749,    351,    681,    842,   1064,
   1674,   3512,   3272,   4474,    686,    603,    232,    919,
    372,    550,    502,    610,    141,    272,    467,    459,
   3072,   6326,   3252,  10678,   1354,   1451,    432,   2016,
   1747,   2537,   1269,   2667,    267,    424,    602,    625,
      0,    841,      0,   1251,   1133,   1275,    342,   1970,
      0,    534,    107,    865,    181,    228,    272,    401},
/* ACA */
{  4300,   2579,   3195,   2735,  70464,  35853,  18993,  36681,
   1840,   4516,   1039,   3463,   3103,   1946,   5607,   2147,
   1277,    520,   2129,    547,   3628,   1634,    854,   2310,
    276,    300,    383,    229,    569,    568,   1783,    601,
   1993,    673,   1966,    983,   9590,   5142,   2048,   5778,
   1320,    923,    844,    726,   2163,   1576,   4056,   1895,
      0,    156,      0,    206,   7353,   4581,   1837,   6050,
      0,    400,    152,    526,    798,    329,    964,    429},
/* ACC */
{  1732,   4969,   1797,   2939,  35853,  72268,  15062,  34881,
    700,   8759,    604,   3844,   1277,   4040,   2951,   2200,
    631,    812,   1389,    497,   1543,   2597,    560,   1934,
    187,    386,    266,    188,    245,    898,   1240,    508,
    960,   1167,   1203,    802,   4073,  10379,   1650,   5231,
    772,   1616,    606,    756,    982,   2920,   2745,   1760,
      0,    298,      0,    241,   3687,   7244,   1275,   5121,
      0,    702,     97,    486,    329,    582,    559,    394},
/* ACG */
{  1129,    943,   1538,    939,  18993,  15062,   9680,  13523,
    428,   1896,    506,   1292,    626,    757,   2625,    774,
    345,    184,    912,    160,    704,    598,    346,    620,
    114,    151,    190,     75,    154,    256,    787,    244,
    495,    296,    748,    285,   2165,   1989,   1035,   1818,
    372,    384,    366,    238,    558,    623,   1746,    707,
      0,     78,      0,     91,   1860,   1833,    795,   1914,
      0,    152,     75,    196,    207,    136,    397,    169},
/* ACT */
{  1860,   3302,   1577,   4546,  36681,  34881,  13523,  50168,
    703,   5687,    538,   5850,   1335,   2174,   3070,   3693,
    743,    527,   1270,    660,   1826,   1601,    581,   2933,
    152,    270,    225,    260,    267,    532,   1170,    758,
   1033,    889,   1236,   1262,   4658,   5276,   1507,   8364,
    812,   1026,    650,   1014,   1034,   1677,   2847,   2631,
      0,    232,      0,    311,   4237,   4379,   1355,   7128,
      0,    498,     74,    697,    389,    391,    616,    571},
/* AGA */
{ 20646,   1338,  15336,   1376,   1840,    700,    428,    703,
  64584,   2232,  35379,   1962,    573,    252,    945,    270,
   2600,   1083,   3685,    966,    402,    236,    173,    267,
  15497,  11485,  15504,   8562,    191,    208,    573,    224,
   1648,    399,   1144,    455,    712,    450,    198,    463,
   2666,    844,   1151,    598,    307,    176,    428,    207,
      0,    145,      0,    141,    429,    331,    121,    422,
      0,    282,    472,    326,    179,     71,    211,    109},
/* AGC */
{  2819,  19096,   2826,  11367,   4516,   8759,   1896,   5687,
   2232, 115828,   2078,  54263,    429,   1496,   1021,    954,
   1079,   2354,   2379,   1395,    793,   1227,    331,    950,
    551,   1609,    766,    673,    130,    512,    632,    353,
   1815,   3818,   2273,   2615,   2075,   4419,    844,   2847,
   4073,  11112,   3476,   4472,    387,   1015,    920,    802,
      0,    637,      0,    450,   2489,   4806,   1038,   3883,
      0,   3222,    187,   2186,    179,    398,    264,    328},
/* AGG */
{ 12905,   1159,  19246,   1093,   1039,    604,    506,    538,
  35379,   2078,  50950,   1647,    317,    220,   1264,    223,
   1629,   1036,   4284,    760,    336,    210,    181,    222,
  10455,  11329,  16771,   6741,    133,    219,    629,    169,
    933,    322,   1571,    353,    516,    430,    228,    344,
   1226,    734,   2008,    458,    173,    190,    577,    164,
      0,    123,      0,    139,    265,    262,    113,    300,
      0,    282,    928,    256,    134,     58,    217,     93},
/* AGT */
{  2447,   9374,   2259,  13749,   3463,   3844,   1292,   5850,
   1962,  54263,   1647,  65410,    369,    634,    792,   1164,
    858,   1161,   1612,   1511,    674,    552,    217,    952,
    382,    611,    555,    753,    105,    230,    416,    332,
   1539,   2036,   1794,   3160,   1780,   2003,    574,   2839,
   3018,   4831,   2430,   5223,    288,    490,    717,    807,
      0,    353,      0,    544,   2045,   2280,    693,   3451,
      0,   1656,    154,   2286,    134,    183,    235,    338},
/* ATA */
{  1031,    309,    649,    351,   3103,   1277,    626,   1335,
    573,    429,    317,    369,  30638,  28175,   7711,  25576,
    288,    130,    452,    155,    462,    224,    121,    301,
    108,     87,    109,     75,   2138,   2486,   5627,   2186,
    523,    109,    445,    168,   1415,    765,    360,    898,
    314,    178,    184,    119,   5479,   4631,   9780,   4914,
      0,    109,      0,    134,    571,    331,    130,    345,
      0,    106,     76,    119,   2544,    701,   2496,    891},
/* ATC */
{   476,   1115,    503,    681,   1946,   4040,    757,   2174,
    252,   1496,    220,    634,  28175, 122704,   7704,  65193,
    264,    333,    536,    242,    339,    490,    134,    372,
     64,    141,    105,     83,   1923,   8057,   8591,   4462,
    335,    304,    395,    174,   1009,   2396,    467,   1451,
    213,    488,    230,    194,   4773,  18081,  15695,  10382,
      0,    402,      0,    282,    392,    642,    129,    507,
      0,    369,     99,    269,   1904,   2777,   2626,   1832},
/* ATG */
{  1993,    772,   3109,    842,   5607,   2951,   2625,   3070,
    945,   1021,   1264,    792,   7711,   7704, 214848,   7959,
    927,    401,   2499,    454,    907,    528,    366,    670,
    227,    290,    444,    189,   3664,   5569,  20280,   4883,
    875,    299,   1505,    324,   2756,   2046,   1059,   2026,
    618,    451,    736,    327,   3533,   3494,  14027,   3874,
      0,    332,      0,    338,   1381,    920,    505,   1193,
      0,    251,    429,    279,   3689,   1481,   8984,   1728},
/* ATT */
{   580,    663,    515,   1064,   2147,   2200,    774,   3693,
    270,    954,    223,   1164,  25576,  65193,   7959,  92314,
    253,    255,    549,    359,    406,    331,    111,    591,
     93,     95,     99,    102,   1876,   4633,   7578,   5795,
    370,    209,    482,    404,   1252,   1323,    412,   2096,
    269,    309,    196,    306,   4791,   9574,  13874,  15034,
      0,    262,      0,    369,    431,    495,    137,    731,
      0,    242,    105,    365,   1920,   1659,   2852,   2625},
/* CAA */
{  7019,   1484,   5156,   1674,   1277,    631,    345,    743,
   2600,   1079,   1629,    858,    288,    264,    927,    253,
  59834,   5592,  80759,   5231,   4070,   1546,   1134,   1926,
   2527,   1646,   2395,   1298,    920,    757,   2285,    834,
   6770,   1219,   5911,   1431,   1215,    741,    369,    798,
    892,    441,    544,    311,    373,    281,    659,    354,
      0,    567,      0,    580,   1137,    755,    264,    917,
      0,    304,    312,    320,    464,    237,    669,    285},
/* CAC */
{  1702,   5608,   1891,   3512,    520,    812,    184,    527,
   1083,   2354,   1036,   1161,    130,    333,    401,    255,
   5592,  99998,  12592,  54421,   1206,   1967,    553,   1426,
   1135,   4128,   1592,   1884,    305,   1413,   1218,    901,
   1120,   2009,   1556,   1296,    426,    753,    185,    486,
    368,    827,    318,    334,    139,    325,    373,    210,
      0,   6304,      0,   3889,    519,    988,    190,    845,
      0,    993,    208,    720,    217,   1193,    348,    849},
/* CAG */
{ 10180,   3282,  15687,   3272,   2129,   1389,    912,   1270,
   3685,   2379,   4284,   1612,    452,    536,   2499,    549,
  80759,  12592, 295886,  10607,   5905,   3254,   2960,   3395,
   3393,   3962,   8200,   2514,   1401,   1989,   7001,   1805,
  10336,   2693,  18331,   2695,   1920,   1873,    919,   1588,
   1203,   1053,   1424,    675,    568,    600,   1676,    652,
      0,   1225,      0,   1035,   1827,   1581,    678,   1639,
      0,    615,   1190,    619,    885,    482,   1801,    496},
/* CAT */
{  1597,   3172,   1526,   4474,    547,    497,    160,    660,
    966,   1395,    760,   1511,    155,    242,    454,    359,
   5231,  54421,  10607,  65564,   1113,   1200,    454,   1972,
    954,   2015,   1215,   2269,    270,    910,   1126,   1287,
   1256,   1251,   1428,   1844,    448,    393,    166,    580,
    367,    494,    278,    508,    127,    205,    338,    261,
      0,   3419,      0,   4605,    492,    614,    151,   1025,
      0,    617,    251,   1012,    236,    684,    352,    954},
/* CCA */
{  1178,    567,    945,    686,   3628,   1543,    704,   1826,
    402,    793,    336,    674,    462,    339,    907,    406,
   4070,   1206,   5905,   1113,  96410,  44179,  22483,  58526,
    465,    384,    509,    283,   1421,    934,   3255,   1141,
   1194,    414,   1309,    514,   4498,   2561,   1031,   2962,
    672,    545,    470,    381,    683,    519,   1253,    566,
      0,    139,      0,    152,   6667,   3396,   1365,   4480,
      0,    155,    137,    230,    741,    245,   1030,    300},
/* CCC */
{   530,    801,    524,    603,   1634,   2597,    598,   1601,
    236,   1227,    210,    552,    224,    490,    528,    331,
   1546,   1967,   3254,   1200,  44179,  76760,  17395,  50743,
    238,    651,    340,    293,    494,   1849,   1646,   1030,
    599,    572,    843,    409,   2153,   4232,    902,   2634,
    412,    684,    384,    325,    307,    695,    761,    451,
      0,    298,      0,    175,   2904,   6600,   1004,   4368,
      0,    290,     55,    215,    292,    503,    488,    308},
/* CCG */
{   358,    232,    430,    232,    854,    560,    346,    581,
    173,    331,    181,    217,    121,    134,    366,    111,
   1134,    553,   2960,    454,  22483,  17395,  12070,  18860,
    162,    224,    371,    130,    347,    378,   1695,    416,
    351,    191,    513,    189,   1157,   1109,    688,    970,
    223,    262,    242,    131,    158,    202,    498,    201,
      0,     75,      0,     61,   1485,   1287,    642,   1472,
      0,     85,     94,     97,    203,     98,    442,    122},
/* CCT */
{   689,    663,    604,    919,   2310,   1934,    620,   2933,
    267,    950,    222,    952,    301,    372,    670,    591,
   1926,   1426,   3395,   1972,  58526,  50743,  18860,  93780,
    273,    413,    312,    541,    576,   1280,   1943,   2255,
    744,    431,    873,    596,   3013,   3191,    934,   4810,
    493,    553,    374,    541,    358,    492,    914,    888,
      0,    219,      0,    305,   3913,   4322,   1159,   8477,
      0,    284,     70,    475,    371,    312,    691,    656},
/* CGA */
{  4714,    333,   4333,    372,    276,    187,    114,    152,
  15497,    551,  10455,    382,    108,     64,    227,     93,
   2527,   1135,   3393,    954,    465,    238,    162,    273,
  15332,  12217,  13400,   9093,    207,    219,    485,    175,
    414,    115,    402,    127,    189,    148,     65,    109,
    488,    210,    307,    184,     73,     57,    129,     57,
      0,    107,      0,    145,    177,    137,     56,    136,
      0,    275,    379,    268,     78,     49,    133,     60},
/* CGC */
{  4108,    942,   5231,    550,    300,    386,    151,    270,
  11485,   1609,  11329,    611,     87,    141,    290,     95,
   1646,   4128,   3962,   2015,    384,    651,    224,    413,
  12217,  40282,  16178,  15122,    124,    750,    658,    364,
    303,    314,    554,    219,    174,    409,    112,    191,
    316,    804,    303,    252,     53,    146,    202,    100,
      0,    439,      0,    280,    153,    397,     85,    239,
      0,   1232,    354,    522,     95,    198,    143,     87},
/* CGG */
{  5511,    541,   7539,    502,    383,    266,    190,    225,
  15504,    766,  16771,    555,    109,    105,    444,     99,
   2395,   1592,   8200,   1215,    509,    340,    371,    312,
  13400,  16178,  26516,  10668,    210,    321,   1172,    299,
    440,    197,    905,    200,    275,    238,    165,    238,
    423,    398,    697,    190,     85,    116,    306,     99,
      0,    149,      0,    141,    215,    238,    114,    214,
      0,    409,   1065,    335,    127,     84,    285,     88},
/* CGT */
{  2693,    433,   2968,    610,    229,    188,     75,    260,
   8562,    673,   6741,    753,     75,     83,    189,    102,
   1298,   1884,   2514,   2269,    283,    293,    130,    541,
   9093,  15122,  10668,  14346,     98,    283,    392,    397,
    288,    147,    313,    212,    164,    190,     63,    213,
    260,    300,    190,    319,     48,     76,    133,    100,
      0,    234,      0,    349,    143,    191,     50,    311,
      0,    474,    218,    748,     73,     96,    132,    138},
/* CTA */
{   363,    126,    281,    141,    569,    245,    154,    267,
    191,    130,    133,    105,   2138,   1923,   3664,   1876,
    920,    305,   1401,    270,   1421,    494,    347,    576,
    207,    124,    210,     98,  14232,  14222,  38640,  13070,
    259,     60,    252,     71,    563,    319,    142,    352,
    134,     72,     87,     59,   1125,    911,   2357,    902,
      0,    183,      0,    206,    621,    278,    149,    336,
      0,     84,    184,    129,  10238,   1246,  12502,   1349},
/* CTC */
{   366,    446,    370,    272,    568,    898,    256,    532,
    208,    512,    219,    230,   2486,   8057,   5569,   4633,
    757,   1413,   1989,    910,    934,   1849,    378,   1280,
    219,    750,    321,    283,  14222,  77258,  64691,  37051,
    269,    220,    393,    139,    561,   1136,    264,    687,
    167,    268,    148,    108,   1260,   4169,   3973,   2309,
      0,    840,      0,    512,    440,    859,    144,    563,
      0,    540,    235,    339,  10893,   6333,  18007,   3656},
/* CTG */
{   884,    471,   1292,    467,   1783,   1240,    787,   1170,
    573,    632,    629,    416,   5627,   8591,  20280,   7578,
   2285,   1218,   7001,   1126,   3255,   1646,   1695,   1943,
    485,    658,   1172,    392,  38640,  64691, 236480,  56553,
    712,    261,   1205,    260,   1777,   1525,    899,   1501,
    353,    352,    498,    231,   3427,   4149,  13822,   4091,
      0,    819,      0,    839,   1564,   1083,    619,   1233,
      0,    439,   1055,    504,  30609,   5669,  65624,   6145},
/* CTT */
{   398,    263,    364,    459,    601,    508,    244,    758,
    224,    353,    169,    332,   2186,   4462,   4883,   5795,
    834,    901,   1805,   1287,   1141,   1030,    416,   2255,
    175,    364,    299,    397,  13070,  37051,  56553,  56660,
    285,    167,    347,    221,    646,    712,    242,    979,
    172,    186,    138,    169,   1209,   2079,   3407,   3272,
      0,    520,      0,    649,    516,    495,    145,    928,
      0,    316,    225,    438,  10125,   3229,  17355,   5183},
/* GAA */
{  8048,   2525,   5401,   3072,   1993,    960,    495,   1033,
   1648,   1815,    933,   1539,    523,    335,    875,    370,
   6770,   1120,  10336,   1256,   1194,    599,    351,    744,
    414,    303,    440,    288,    259,    269,    712,    285,
 206378,  24935, 169440,  28911,   4616,   2196,   1162,   2462,
   5785,   1997,   2920,   1599,   1082,    645,   1894,    798,
      0,    234,      0,    280,   1146,    745,    302,    981,
      0,    137,    126,    187,    275,    130,    423,    180},
/* GAC */
{  1636,  10085,   1605,   6326,    673,   1167,    296,    889,
    399,   3818,    322,   2036,    109,    304,    299,    209,
   1219,   2009,   2693,   1251,    414,    572,    191,    431,
    115,    314,    197,    147,     60,    220,    261,    167,
  24935, 179326,  32603, 109888,   1258,   2641,    589,   1529,
   2157,   4864,   1555,   2004,    237,    701,    547,    465,
      0,    608,      0,    406,    570,    971,    216,    755,
      0,    263,     45,    221,     87,    175,    116,    114},
/* GAG */
{  5765,   3034,   9352,   3252,   1966,   1203,    748,   1236,
   1144,   2273,   1571,   1794,    445,    395,   1505,    482,
   5911,   1556,  18331,   1428,   1309,    843,    513,    873,
    402,    554,    905,    313,    252,    393,   1205,    347,
 169440,  32603, 306080,  32897,   4089,   2938,   2087,   2798,
   3762,   2387,   4836,   1598,    976,    823,   3117,    849,
      0,    290,      0,    328,   1111,    994,    398,   1118,
      0,    214,    240,    202,    279,    167,    507,    182},
/* GAT */
{  1836,   6768,   1732,  10678,    983,    802,    285,   1262,
    455,   2615,    353,   3160,    168,    174,    324,    404,
   1431,   1296,   2695,   1844,    514,    409,    189,    596,
    127,    219,    200,    212,     71,    139,    260,    221,
  28911, 109888,  32897, 182052,   1528,   1598,    520,   2698,
   2256,   2810,   1611,   3626,    295,    412,    647,    868,
      0,    379,      0,    709,    680,    670,    209,   1178,
      0,    186,     54,    334,     86,     98,    156,    206},
/* GCA */
{  1979,   1131,   1492,   1354,   9590,   4073,   2165,   4658,
    712,   2075,    516,   1780,   1415,   1009,   2756,   1252,
   1215,    426,   1920,    448,   4498,   2153,   1157,   3013,
    189,    174,    275,    164,    563,    561,   1777,    646,
   4616,   1258,   4089,   1528,  71446,  42285,  15562,  48702,
   4400,   2448,   2627,   1992,   3498,   2533,   7151,   3146,
      0,    171,      0,    203,   7531,   4721,   1802,   6059,
      0,    480,    146,    584,    774,    328,    948,    484},
/* GCC */
{  1019,   2290,   1138,   1451,   5142,  10379,   1989,   5276,
    450,   4419,    430,   2003,    765,   2396,   2046,   1323,
    741,    753,   1873,    393,   2561,   4232,   1109,   3191,
    148,    409,    238,    190,    319,   1136,   1525,    712,
   2196,   2641,   2938,   1598,  42285, 109722,  16086,  57041,
   2770,   5134,   2293,   2068,   1890,   6486,   5502,   3573,
      0,    351,      0,    255,   4869,   9986,   1800,   6924,
      0,   1065,    147,    697,    421,    792,    700,    509},
/* GCG */
{   494,    489,    599,    432,   2048,   1650,   1035,   1507,
    198,    844,    228,    574,    360,    467,   1059,    412,
    369,    185,    919,    166,   1031,    902,    688,    934,
     65,    112,    165,     63,    142,    264,    899,    242,
   1162,    589,   2087,    520,  15562,  16086,   9674,  15220,
   1135,   1096,   1337,    655,    841,   1072,   3562,   1117,
      0,    103,      0,     73,   1755,   1851,    959,   2061,
      0,    229,     72,    215,    164,    169,    447,    166},
/* GCT */
{  1188,   1536,   1064,   2016,   5778,   5231,   1818,   8364,
    463,   2847,    344,   2839,    898,   1451,   2026,   2096,
    798,    486,   1588,    580,   2962,   2634,    970,   4810,
    109,    191,    238,    213,    352,    687,   1501,    979,
   2462,   1529,   2798,   2698,  48702,  57041,  15220,  87906,
   2987,   2955,   2050,   3068,   2160,   3308,   5558,   5958,
      0,    277,      0,    360,   5339,   5936,   1674,   9628,
      0,    697,     92,    943,    476,    502,    774,    796},
/* GGA */
{  1686,   1728,   1215,   1747,   1320,    772,    372,    812,
   2666,   4073,   1226,   3018,    314,    213,    618,    269,
    892,    368,   1203,    367,    672,    412,    223,    493,
    488,    316,    423,    260,    134,    167,    353,    172,
   5785,   2157,   3762,   2256,   4400,   2770,   1135,   2987,
 107006,  49273,  52181,  36526,    861,    482,   1260,    577,
      0,     73,      0,    116,   1063,    876,    290,   1030,
      0,    415,    258,    494,    149,    113,    209,    168},
/* GGC */
{   799,   4186,    885,   2537,    923,   1616,    384,   1026,
    844,  11112,    734,   4831,    178,    488,    451,    309,
    441,    827,   1053,    494,    545,    684,    262,    553,
    210,    804,    398,    300,     72,    268,    352,    186,
   1997,   4864,   2387,   2810,   2448,   5134,   1096,   2955,
  49273,  97496,  36303,  38034,    351,   1106,    910,    657,
      0,    251,      0,    162,    815,   1434,    371,   1134,
      0,   1362,    192,    854,     90,    226,    180,    184},
/* GGG */
{   938,   1544,   1386,   1269,    844,    606,    366,    650,
   1151,   3476,   2008,   2430,    184,    230,    736,    196,
    544,    318,   1424,    278,    470,    384,    242,    374,
    307,    303,    697,    190,     87,    148,    498,    138,
   2920,   1555,   4836,   1611,   2627,   2293,   1337,   2050,
  52181,  36303,  57212,  24188,    525,    471,   1707,    502,
      0,     81,      0,     97,    758,    688,    323,    720,
      0,    381,    636,    380,    116,    105,    291,    110},
/* GGT */
{   654,   1992,    628,   2667,    726,    756,    238,   1014,
    598,   4472,    458,   5223,    119,    194,    327,    306,
    311,    334,    675,    508,    381,    325,    131,    541,
    184,    252,    190,    319,     59,    108,    231,    169,
   1599,   2004,   1598,   3626,   1992,   2068,    655,   3068,
  36526,  38034,  24188,  45116,    310,    419,    741,    822,
      0,    105,      0,    183,    609,    662,    219,   1035,
      0,    566,    110,    917,     83,     95,    117,    171},
/* GTA */
{   591,    242,    432,    267,   2163,    982,    558,   1034,
    307,    387,    173,    288,   5479,   4773,   3533,   4791,
    373,    139,    568,    127,    683,    307,    158,    358,
     73,     53,     85,     48,   1125,   1260,   3427,   1209,
   1082,    237,    976,    295,   3498,   1890,    841,   2160,
    861,    351,    525,    310,  18232,  13787,  36807,  15198,
      0,     85,      0,     96,    728,    403,    163,    534,
      0,    131,     64,    157,   1632,    440,   1673,    557},
/* GTC */
{   364,    636,    370,    424,   1576,   2920,    623,   1677,
    176,   1015,    190,    490,   4631,  18081,   3494,   9574,
    281,    325,    600,    205,    519,    695,    202,    492,
     57,    146,    116,     76,    911,   4169,   4149,   2079,
    645,    701,    823,    412,   2533,   6486,   1072,   3308,
    482,   1106,    471,    419,  13787,  54440,  44453,  27897,
      0,    332,      0,    208,    589,    988,    175,    715,
      0,    377,     77,    272,    972,   1927,   1567,   1115},
/* GTG */
{   976,    590,   1230,    602,   4056,   2745,   1746,   2847,
    428,    920,    577,    717,   9780,  15695,  14027,  13874,
    659,    373,   1676,    338,   1253,    761,    498,    914,
    129,    202,    306,    133,   2357,   3973,  13822,   3407,
   1894,    547,   3117,    647,   7151,   5502,   3562,   5558,
   1260,    910,   1707,    741,  36807,  44453, 162686,  46454,
      0,    303,      0,    278,   1444,   1089,    539,   1354,
      0,    391,    294,    478,   2622,   1431,   5885,   1606},
/* GTT */
{   468,    466,    389,    625,   1895,   1760,    707,   2631,
    207,    802,    164,    807,   4914,  10382,   3874,  15034,
    354,    210,    652,    261,    566,    451,    201,    888,
     57,    100,     99,    100,    902,   2309,   4091,   3272,
    798,    465,    849,    868,   3146,   3573,   1117,   5958,
    577,    657,    502,    822,  15198,  27897,  46454,  48518,
      0,    181,      0,    279,    726,    727,    192,   1127,
      0,    318,    112,    377,   1102,   1049,   1825,   1770},
/* TAA */
{     0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
#ifdef INCLUDE_STOP_CODON_FREQUENCIES
   2564,      0,    732,      0,      0,      0,      0,      0,
   1164,      0,      0,      0,      0,      0,      0,      0},
#else /* INCLUDE_STOP_CODON_FREQUENCIES */
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0},
#endif /* INCLUDE_STOP_CODON_FREQUENCIES */
/* TAC */
{   247,   1308,    273,    841,    156,    298,     78,    232,
    145,    637,    123,    353,    109,    402,    332,    262,
    567,   6304,   1225,   3419,    139,    298,     75,    219,
    107,    439,    149,    234,    183,    840,    819,    520,
    234,    608,    290,    379,    171,    351,    103,    277,
     73,    251,     81,    105,     85,    332,    303,    181,
      0, 138828,      0,  75342,    368,   1128,    149,    794,
      0,   2235,    555,   1461,    349,  11866,    429,   8131},
/* TAG */
{     0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
#ifdef INCLUDE_STOP_CODON_FREQUENCIES
    732,      0,   1132,      0,      0,      0,      0,      0,
    534,      0,      0,      0,      0,      0,      0,      0},
#else /* INCLUDE_STOP_CODON_FREQUENCIES */
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0},
#endif /* INCLUDE_STOP_CODON_FREQUENCIES */
/* TAT */
{   261,    803,    257,   1251,    206,    241,     91,    311,
    141,    450,    139,    544,    134,    282,    338,    369,
    580,   3889,   1035,   4605,    152,    175,     61,    305,
    145,    280,    141,    349,    206,    512,    839,    649,
    280,    406,    328,    709,    203,    255,     73,    360,
    116,    162,     97,    183,     96,    208,    278,    279,
      0,  75342,      0,  95904,    406,    708,    137,   1247,
      0,   1327,    535,   2349,    352,   7329,    417,  10375},
/* TCA */
{  1159,   1077,    919,   1133,   7353,   3687,   1860,   4237,
    429,   2489,    265,   2045,    571,    392,   1381,    431,
   1137,    519,   1827,    492,   6667,   2904,   1485,   3913,
    177,    153,    215,    143,    621,    440,   1564,    516,
   1146,    570,   1111,    680,   7531,   4869,   1755,   5339,
   1063,    815,    758,    609,    728,    589,   1444,    726,
      0,    368,      0,    406,  53222,  30645,  12204,  36223,
      0,    868,    374,   1086,   1891,    603,   1941,    840},
/* TCC */
{   733,   2021,    702,   1275,   4581,   7244,   1833,   4379,
    331,   4806,    262,   2280,    331,    642,    920,    495,
    755,    988,   1581,    614,   3396,   6600,   1287,   4322,
    137,    397,    238,    191,    278,    859,   1083,    495,
    745,    971,    994,    670,   4721,   9986,   1851,   5936,
    876,   1434,    688,    662,    403,    988,   1089,    727,
      0,   1128,      0,    708,  30645,  79396,  11832,  45278,
      0,   2263,    247,   1381,    655,   1997,   1030,   1070},
/* TCG */
{   292,    360,    298,    342,   1837,   1275,    795,   1355,
    121,   1038,    113,    693,    130,    129,    505,    137,
    264,    190,    678,    151,   1365,   1004,    642,   1159,
     56,     85,    114,     50,    149,    144,    619,    145,
    302,    216,    398,    209,   1802,   1800,    959,   1674,
    290,    371,    323,    219,    163,    175,    539,    192,
      0,    149,      0,    137,  12204,  11832,   6594,  12964,
      0,    325,    220,    343,    275,    202,    772,    285},
/* TCT */
{   981,   1682,    830,   1970,   6050,   5121,   1914,   7128,
    422,   3883,    300,   3451,    345,    507,   1193,    731,
    917,    845,   1639,   1025,   4480,   4368,   1472,   8477,
    136,    239,    214,    311,    336,    563,   1233,    928,
    981,    755,   1118,   1178,   6059,   6924,   2061,   9628,
   1030,   1134,    720,   1035,    534,    715,   1354,   1127,
      0,    794,      0,   1247,  36223,  45278,  12964,  78934,
      0,   1651,    257,   2645,    849,   1209,   1299,   2209},
/* TGA */
{     0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0,
#ifdef INCLUDE_STOP_CODON_FREQUENCIES
   1164,      0,    534,      0,      0,      0,      0,      0,
   3580,      0,      0,      0,      0,      0,      0,      0},
#else /* INCLUDE_STOP_CODON_FREQUENCIES */
      0,      0,      0,      0,      0,      0,      0,      0,
      0,      0,      0,      0,      0,      0,      0,      0},
#endif /* INCLUDE_STOP_CODON_FREQUENCIES */
/* TGC */
{   280,    817,    311,    534,    400,    702,    152,    498,
    282,   3222,    282,   1656,    106,    369,    251,    242,
    304,    993,    615,    617,    155,    290,     85,    284,
    275,   1232,    409,    474,     84,    540,    439,    316,
    137,    263,    214,    186,    480,   1065,    229,    697,
    415,   1362,    381,    566,    131,    377,    391,    318,
      0,   2235,      0,   1327,    868,   2263,    325,   1651,
      0, 103986,    771,  57765,    166,   1430,    279,   1086},
/* TGG */
{   264,    100,    441,    107,    152,     97,     75,     74,
    472,    187,    928,    154,     76,     99,    429,    105,
    312,    208,   1190,    251,    137,     55,     94,     70,
    379,    354,   1065,    218,    184,    235,   1055,    225,
    126,     45,    240,     54,    146,    147,     72,     92,
    258,    192,    636,    110,     64,     77,    294,    112,
      0,    555,      0,    535,    374,    247,    220,    257,
      0,    771, 189124,    755,    243,    629,   1053,    730},
/* TGT */
{   298,    575,    290,    865,    526,    486,    196,    697,
    326,   2186,    256,   2286,    119,    269,    279,    365,
    320,    720,    619,   1012,    230,    215,     97,    475,
    268,    522,    335,    748,    129,    339,    504,    438,
    187,    221,    202,    334,    584,    697,    215,    943,
    494,    854,    380,    917,    157,    272,    478,    377,
      0,   1461,      0,   2349,   1086,   1381,    343,   2645,
      0,  57765,    755,  88088,    203,    944,    334,   1600},
/* TTA */
{   341,    139,    263,    181,    798,    329,    207,    389,
    179,    179,    134,    134,   2544,   1904,   3689,   1920,
    464,    217,    885,    236,    741,    292,    203,    371,
     78,     95,    127,     73,  10238,  10893,  30609,  10125,
    275,     87,    279,     86,    774,    421,    164,    476,
    149,     90,    116,     83,   1632,    972,   2622,   1102,
      0,    349,      0,    352,   1891,    655,    275,    849,
      0,    166,    243,    203,  25300,   2266,  18809,   2526},
/* TTC */
{   143,    395,    120,    228,    329,    582,    136,    391,
     71,    398,     58,    183,    701,   2777,   1481,   1659,
    237,   1193,    482,    684,    245,    503,     98,    312,
     49,    198,     84,     96,   1246,   6333,   5669,   3229,
    130,    175,    167,     98,    328,    792,    169,    502,
    113,    226,    105,     95,    440,   1927,   1431,   1049,
      0,  11866,      0,   7329,    603,   1997,    202,   1209,
      0,   1430,    629,    944,   2266, 167022,   3062,  90000},
/* TTG */
{   442,    238,    530,    272,    964,    559,    397,    616,
    211,    264,    217,    235,   2496,   2626,   8984,   2852,
    669,    348,   1801,    352,   1030,    488,    442,    691,
    133,    143,    285,    132,  12502,  18007,  65624,  17355,
    423,    116,    507,    156,    948,    700,    447,    774,
    209,    180,    291,    117,   1673,   1567,   5885,   1825,
      0,    429,      0,    417,   1941,   1030,    772,   1299,
      0,    279,   1053,    334,  18809,   3062,  44874,   3690},
/* TTT */
{   172,    305,    170,    401,    429,    394,    169,    571,
    109,    328,     93,    338,    891,   1832,   1728,   2625,
    285,    849,    496,    954,    300,    308,    122,    656,
     60,     87,     88,    138,   1349,   3656,   6145,   5183,
    180,    114,    182,    206,    484,    509,    166,    796,
    168,    184,    110,    171,    557,   1115,   1606,   1770,
      0,   8131,      0,  10375,    840,   1070,    285,   2209,
      0,   1086,    730,   1600,   2526,  90000,   3690, 156062}
 };


static gint CodonSubmat_round(gdouble num){
    return ( num < 0 ) ? num - 0.5 : num + 0.5;
    }
/* Not perfect, but portable */

CodonSubmat *CodonSubmat_create(){
    register CodonSubmat *cs = g_new(CodonSubmat, 1);
    register guchar *alphabet = (guchar*)"ACGTN";
    register gint i, j, a, b, c, freq, m, n, total = 0;
    register gdouble ratio, logbase = log(10.0), prob, chance;
    gint codon_freq_acgtn[CODON_ALPHABETSIZE_ACGTN] = {0};
    gint red_codon[CODON_ALPHABETSIZE_ACGT][8];
    gint count_acgtn[CODON_ALPHABETSIZE_ACGTN]
                    [CODON_ALPHABETSIZE_ACGTN] = {{0}};
    cs->ref_count = 1;
    /* Fill the base index */
    for(i = 0; i < ALPHABETSIZE; i++)
        cs->base_index[i] = 4;
    for(i = 0; i < 4; i++)
        cs->base_index[alphabet[i]]
      = cs->base_index[tolower(alphabet[i])] = i;
    /* Fill the codon index */
    i = 0;
    for(a = 0; a < 5; a++)
        for(b = 0; b < 5; b++)
            for(c = 0; c < 5; c++)
                cs->codon_index[a][b][c] = i++;
    /* Fill the redundant codon table
     * eg. [AAA] -> {AAN,ANA,NAA,ANN,NAN,NNA,NNN}
     */
    i = 0;
    for(a = 0; a < 4; a++)
        for(b = 0; b < 4; b++)
            for(c = 0; c < 4; c++){
                red_codon[i][0] = cs->codon_index[a][b][c];
                red_codon[i][1] = cs->codon_index[a][b][4];
                red_codon[i][2] = cs->codon_index[a][4][c];
                red_codon[i][3] = cs->codon_index[4][b][c];
                red_codon[i][4] = cs->codon_index[a][4][4];
                red_codon[i][5] = cs->codon_index[4][b][4];
                red_codon[i][6] = cs->codon_index[4][4][c];
                red_codon[i][7] = cs->codon_index[4][4][4];
                i++;
                }
    /* Calculate codon frequencies */
    for(i = 0; i < CODON_ALPHABETSIZE_ACGT; i++){
        for(j = 0; j < CODON_ALPHABETSIZE_ACGT; j++){
            codon_freq_acgtn[red_codon[i][0]]
                += local_codon_substitution_count[i][j];
            total += local_codon_substitution_count[i][j];
            }
        }
    /* Calculate redundant codon frequencies */
    for(i = 0; i < CODON_ALPHABETSIZE_ACGT; i++)
        for(j = 1; j < 8; j++)
            codon_freq_acgtn[red_codon[i][j]]
                += codon_freq_acgtn[red_codon[i][0]];
    /* Fill the combined count matrix including redundancies */
    for(i = 0; i < CODON_ALPHABETSIZE_ACGT; i++){
        for(j = 0; j < CODON_ALPHABETSIZE_ACGT; j++){
            freq = local_codon_substitution_count[i][j];
            for(m = 0; m < 8; m++)
                for(n = 0; n < 8; n++)
                    count_acgtn[red_codon[i][m]]
                               [red_codon[j][n]] += freq;
            }
        }
    /* Calculate log odds scores */
    for(i = 0; i < CODON_ALPHABETSIZE_ACGTN; i++){
        for(j = 0; j < CODON_ALPHABETSIZE_ACGTN; j++){
            freq = count_acgtn[i][j];
            if(freq){
                prob =  (gdouble)(freq << 1)
                     / ((gdouble) codon_freq_acgtn[i]
                                + codon_freq_acgtn[j]);
                chance = ((gdouble) codon_freq_acgtn[i]
                                  + codon_freq_acgtn[j])
                       / (gdouble)(total << 1);
                ratio = prob / chance;
                cs->codon_submat[i][j]
                    = CodonSubmat_round(10.0 * (log(ratio)/logbase));
            } else {
                cs->codon_submat[i][j] = -50;
                }
            }
        }
    return cs;
    }

void CodonSubmat_destroy(CodonSubmat *cs){
    if(--cs->ref_count)
        return;
    g_free(cs);
    return;
    }

CodonSubmat *CodonSubmat_share(CodonSubmat *cs){
    cs->ref_count++;
    return cs;
    }

gint CodonSubmat_max_score(CodonSubmat *cs){
    register gint i, j, max;
    g_assert(cs);
    max = cs->codon_submat[0][0];
    for(i = 0; i < CODON_ALPHABETSIZE_ACGTN; i++)
        for(j = 0; j < CODON_ALPHABETSIZE_ACGTN; j++)
            if(max < cs->codon_submat[i][j])
                max = cs->codon_submat[i][j];
    return max;
    }

void CodonSubmat_add_nucleic(CodonSubmat *cs, Submat *nucleic){
    register gint a, b, c, d, e, f, codon_a, codon_b;
    register gchar *alphabet = "ACGTN";
    register gint score;
    for(a = 0; a < 5; a++){
        for(b = 0; b < 5; b++){
            for(c = 0; c < 5; c++){
                for(d = 0; d < 5; d++){
                    for(e = 0; e < 5; e++){
                        for(f = 0; f < 5; f++){
                            score = Submat_lookup(nucleic, alphabet[a],
                                                           alphabet[d])
                                  + Submat_lookup(nucleic, alphabet[b],
                                                           alphabet[e])
                                  + Submat_lookup(nucleic, alphabet[c],
                                                           alphabet[f]);
                            codon_a = cs->codon_index[a][b][c];
                            codon_b = cs->codon_index[d][e][f];
                            cs->codon_submat[codon_a][codon_b] /= 2;
                            cs->codon_submat[codon_a][codon_b] += score;
                            }
                        }
                    }
                }
            }
        }
    return;
    }

